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en : SUMMARY 

o 

(N 



The understanding of electrokinetics for unsaturated conditions is crucial for numerous of 
- r-^ _ geophysical data interpretation. Nevertheless, the behaviour of the streaming potential 

c^ ! coefficient C as a function of the water saturation 5*^ is still discussed. We propose 

here to model both the Richards equation for hydrodynamics and the Poisson's equation 
for electrical potential for unsaturated conditions using ID flnite element method. The 
equations are flrst presented and the numerical scheme is then detailed for the Poisson's 
equation. Then, computed streaming potentials (SP) are compared to recently published 
SP measurements carried out during drainage experiment in a sand column. We show 
that the apparent measurement of AV/AP for the dipoles can provide the streaming 
potential coefficient in these conditions. Two tests have been performed using existing 
models for the SP coefficient and a third one using a new relation. The results show 
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that existing models of unsaturated SP coefficients C{Sw) provide poor results in terms 
of SP magnitude and behaviour. We demonstrate that the unsaturated SP coefficient 
can be until one order of magnitude larger than Cgat, its value at saturation. We finally 
prove that the SP coefficient follows a non-monotonous behaviour with respect to water 
saturation. 

1 INTRODUCTION 



Electric and electromagnetic methods are used in a large range of geophysical applications 
because of their sensitivity to fluids within the crust. The electrical resistivity can be re- 



lated to the permeability arid to the deformation, in 



conditions (iDoussan fc Ruy 



2009 



Henry et al. 



2003 



uU-saturated or in partially -saturated 



Jouniaux et al. 



1994 



20061 ). The self- 



potentials have also been successfully used in the last decades. Self-potentials have been 



observed to detect contaminant p 



trochemical effects (JNaudet et al. 



umes or salted fronts throu g 



2003 



Maineult et al 



2006b 



1 the interpretation of elec- 
a|) . Most of the self-potential 



observations are interpreted through the electrokin etic effect. It has be e n proposed to use 



this electrokinetic effect for the seismic prediction 



Fenoglio et al. 



1995 



Pozzi fc Jouniaux 



19941 ). For hydrological ap plications, some hydraulic properties can 
self-potential observat ions (iGibert fc Pessell l2001l : I Sailhac et al 



3e inferred from the 



2004 



Glover et al 



Glover fc Walker 



2006; 



20091 ) . The W-shaped anomalies classically observed on active volcanoes 



are used to characterize geothermal circulatio ns (IFinizola et al 



2004 



Mauri et al. 



2010; 



Onizawa et al. 



2002, 



2004 



Saracco et al 



20091 ). It has also been proposed to use the self- 



potential monitoring to detect at distance the propagation of a water front in a reservoir 



(ISaunders et al. 



20081 ). Moreover self-potentials have been monitor ed during hydraulic tests 



in boreholes leading to some relationship with the microseismicity (JDarnet et al. 



2006|), and 



showing a non-linear behaviour that could be related to the saturation and desaturation pro- 



cesses (IMaineult et al. 
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20081 ). Streaming potential (SP) results from the coupling between 
fluid (water) flow and electrical current, through the motion of ionic charges of water in the 
pore space. The distribution of ions near the matrix surface is described by the electric dou- 
ble layer, including the diffuse layer in which the num ber of counterions exceeds the number 



of cations adsorbed to the matrix (JDavis et al 



19781 ). The streaming current is caused by 



the motion of ions from the diffuse layer, coming from a pressure difference. This current is 
then balanced by a conduction current leading to the SP. 

In steady state flows through homogeneous media, one can define the streaming potential co- 
efficient C as the ratio between the measured electrical potential difference AV and the driv- 



ing pore water pressure difference AP (jOverbeek 



19521 ). The strea ming potential is a func 



tion of various pa r ameters, and its depend ence on water salinity ( jlshido &: Mizutani 



Jaafar et al 



991 



2009 



Lome et al. 



( JTosha et al. 



Vinogradov et al. 



1981 



20101). water electrical conductivity 



19991 ). pH f llshido fc Mizutani 



1981 



Guichet et al 



IPride fc Morgan 



20061 ) or temperature 



20031) is still studied. 

However, one is forced to notice a lake of data concerning the SP coefficient dependence on 
water conten t. For shallow surface geophysical applications, inclu ding also the seismoelectri c 



conversions ( jDupuis et al. 



20071 ) which is studied in laboratory (JBordes et al. 



2006 



20081), 



a better understanding of electrokinetics for unsaturated conditions is needed. A way for 
the understanding of such a phenomena is to study the streaming potential coefficient as a 
function of saturation. 



Guichet et al 



fl2003[ ). 



The first experimental SP coefficient measurements were performed by 

These authors measured SP during drainage experiments performed by injecting inert gas 

through sand, and inferred a linear relation between the relative SP coeffici ent Cr (i.e C nor- 



malize d by its value at saturation Cgat) and effective water saturation S^- 



Perrier &: Morat 



( I2OOOI ) proposed an empirical expression to explain the dependence of Cr on water content 
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based on a relative permeability model. The implicit assumption was that the electrical 

currents ar e affec ted by unsaturated state in a comparable way than hydrological ffow. 



Reviletal 



(120071 1 proposed recently an other formu l a to c haracterize this dependence also 



based on a relative permeability model. 



Lindeetal 



( 120071 ) proposed another expression to 



model some SP me asurements pe r 



tions to those from 



Allegre et al. 



brme d during a drainage experiment, with similar condi- 



(120 lOf ). However, these studies do not provide a combined 



hydrodynamic a nd electrical app r oach to model the data, as it should be done. For reser- 



voir applications 



Saunders et al. 



(120081) proposed a linear expression for Cr in the case of 



oil imbibition. All the se models pre d ict a monotonous decrease of Cr with decreasing water 



saturation. Recently, 



Allegre et al. 



( I2OIOI ) proposed original SP measurements performed 
during a drainage experiment and measured the first continuous recordings of the SP coef- 
ficient as a function of water saturation. They observed that the SP coefficient exhibits two 
different behaviours as the water saturation decreases. Values of Cj. first increase for decreas- 
ing saturation in the range 0.55 — 0.8 < 5^ < 1, and then decrease from S^ = 0.55 — 0.8 
to residual water saturation. This behaviour was never reported before and called for new 
interpretations of electrokinetic phenomena for unsaturated conditions. 
Streaming potentia ls have been successfully modelled for aquifer propert i es determination 



(Darnet et al.l 



2OO3I ) or for water infiltration conditions (jSailhac et al. 



2004f ). 



Sheffer fc Oldenburg 



( I2OO7I ) proposed a 3D m odelling of streaming potential at the field scale for saturated con- 



ditions. iJacksonI ( I2OIOI ) used a bundle capillary model to compute the SP coefficent as a 
function of water-saturation. He showed that the behaviour of the SP coefficient depends 
on the capillary size distribution, the wetting behaviour of the capillaries, and whether we 
invoke the thin or thick electrical double layer assumption. Depending upon the chosen value 
of the saturation exponent and the irreductible water-saturation, the relative SP coefficient 



may increase at partial saturation, before decreasing to zero at the irreductible saturation. 
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Finally, no previous study took into account hydrodynamics and electrical potential equa- 
tions together, even in a simple geometry, to model the streaming potential coefficient for 
unsaturated conditions. Thus, we propose here to model SP by solving both the Richards 
equation for hydrodynamics and the Poisson's equation for electrical potential using ID finite 
element method. Existing models which describe the behaviour of Cr as a func tion of Sp, are 



tested and compared to a new expression inferred from SP measurements by 



AUegre et al 



(120 lOl ). Thus, after introducing governing equations, computed SP using these models are 



presented and compared to measurements. The results lead to the conclusions that 1) a 
non-monotonous behaviour of Cr is required to fit the measurements; 2) The apparent mea- 
surement of AV/AP for the dipoles can provide the streaming potential coefficient in these 
conditions. 



2 GOVERNING EQUATIONS 



2.1 Hydrodynamic equations 



Combining the mass conservation equ ation to t 



le ID generalized Darcy's law leads to the 



mixed form of the Richards equation (JRichardi 

porous media, 

de{h) d 



193l[ ). which describes unsaturated flow in 



A-('»(|-l 







(1) 



dt dz 
where 9{h) is the volumetric water content [m'^.m"^], dependent on the pressure head h [m]. 

The parameter K^ which is also a function of the pressure head, is the hydraulic conductivity 

[m.s~^], t is time [s], and z is the vertical coordinate [m] taken to be positive downward. 

The dependence of the hydraulic conductivity and pressure head on water content is non- 



linear. Numerous retention and relative permeabi 



this dependence (I Gardner 



1958 



Brooks fc Corey 



ity models are able to take into account 



1964 



van Genuchten 



1980). The models 
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which have been chosen for this work were proposed respectively by lBrooks fc Coreyl (119641 ). 

I ( '^'^]^ 

— Or 

hn 



Se 



ha\'^ ^ ha 

\h\ I \h\ 



Uk — Or 



(2) 



and 



MualemI (119761 ). 



1 ,if nJ>l 

\h\ 



K{Se) = Ks.S, 



L+2+i 



(3) 



with 5*6 the effective water saturation, Og the water content at saturation [m^.ni"^], also equal 
to porosity 0, 6r the residual water content [m'^.ni"^], and Kg the hydraulic conductivity 
at saturation [m.s^^]. The effective water saturation can be expressed by: Se = {Suj — 
S'^)/(l — S*^), with Sw {Sw = 9/(j)) and S"^ (5^ = Or/(p) the water saturation and the 
residual water saturation respectively. The parameter A in equation [2] is a measure of the 
pore size distribution and characterizes the medium granulometry. Thus, higher the value 
of A is, higher hom ogeneous the medium i s. The second hydrodynamic parameter h^ is the 



air entry pressure ( iBrooks fc Corey 



1964J ). The last parameter L takes into acc ount the 



tortuosity and is chosen as L = 0.5, which is a common value in the literature (JMualem 



19761). 



Some initial and boundary conditions are necessary to solve the Richards equation. Initial 
condition is a saturated sand at the hydrostatic equilibrium, so that the water pressures 
(expressed in cm of water) at the top and the bottom boundary of the system are respectively 
h{z = 0) = and h{z = I) = I. Considering the drainage experiment which will be presented 
in this work, two conditions have been used. The first one is a zero flux qo{t) at the top 
of the system (Neumann type) and the second one is a constant pressure head hi (t) at the 
bottom (Dirichlet type). These conditions are written as. 



dh 
h{z = l,t) = hi{t) and (-K{h)— + K 



qo{t) 



(4) 



2 = 
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where z = oi z = I and / the length of the syste m. The Richards equa tion has been 



1913), with a fully 



discretized using the Galerkin finite element method ( iPinder &: Gray 

implicit scheme in time. To take into account dependency between h, 6 and K, the equation 

is linearized using the Newton-Raphson method. This approach has bee n used for decades to 



solve t his equation, and the detailed scheme can be found for example in 



Lehmann fc Ackerer 



(11998h 



2.2 Electrokinetic theory 



Coupled fluxes can be described by the general equation, 

N 



(5) 



whi ch link the for ces X,- to macroscopic fluxes Jj, through transport coupling coefficients 



Cij (lOnsager 



193ll ). The global SP field can be described as the sum of several contributions 
creating electrical current sources. These current sources derive from macroscopic potentials, 
through electrochemical effects (e.g. concentration gradients), electrokinetics (e.g. electrical 
potential gradients, electro-osmosis) or thermo-electrical effects (e.g. temperature gradients). 
Considering some of these potential fields, the total electrical current density can be written 



as: 



J = -£ 



VT 



T" 



T 



+ V$j + C,VV - C,kVP, 



(6) 



where T is the temperature [K], Ct = —T-rgt [W.m.s~^] is given by the Peltier effect, 
$j is the junction potential [V], which also can be expressed as a function of chemical 
concentration gradients in el ectrolytes CVCrJCr) by V^j = arn.'VCrJCr., with am. the flui d 



junction coupling coefficient ( jNaudet et al. 



2003 



Maineult et al. 



2008 : 



Jouniaux et al. 



20091). 



The parameter V is the electrical potential [V], thus Ohm's law identify Cg = — cr^, with 0"^ 
the bulk electrical conductivity [S.m^^]. The last term of equation E] describes electrokinetic 
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effects created by the driving pressure gradient V-P , through the electrokinetic couphng £efc, 



defined as Cek = —cTrC [A. Pa .m ] ( iPride 



19941 ) ■ The parameter C [V.Pa ^] is known as 



the SP coefficient. The total water pressure P can be inferred from water pressures h with: 
P = Pwg{h — z), where p^ is the water density [kg.m~^], g is gravity and z is the vertical 
location taken to be positive downward. 

Considering a constant temperature, and no concentration gradients, one can write the 
following coupled equation: 

J = CeW - Cek'^P, (7) 

or, 

J = -arW + arCVP. (8) 

Without any external current sources, the conservation of the total current density implies, 

V ■ J = 0. (9) 

In the case of heterogeneous medium, one can assume for example a tabular medium with 
electrical conductivity and SP coefficient contrasts, then equation [8] through [9] leads to the 
following Poisson's equation: 

V ■ J = -(T,VV - W ■ Var + V(C(T,) • VP + arCV^P = 0. (10) 

In addition to primary current sources occuring from the term arCV P, some secondary 
sources linked to the term 'V{Car) ■ VP appear. These sources are located at boundaries 
formed by electrical conductivity and SP coefficient constrasts. 

In the case of an homogeneous medium and without any contrasts of a,, and C, the equation 
[To] reduces to: 

VV = CV'^P. (11) 
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Considering for example a steady-state saturated flow through a capillary the SP coefficient 

C [V.Pa^^] can be expressed as the ratio between the SP difference AV [V] and the driving- 
pressure difference AP [Pa], 

AV 



C 



AP' 



(12) 



In one dimension, the equation [H] can be written as. 



d_ 

dz 



-(Jr 



dV 
dz 



d_ 

dz 






0. 



(13) 



2.3 Discretization of the Poisson's equation 



The Poisson's equation was solved using ID finite element method. A first order basis func- 
tion (f){z) has been chosen to discretize the Poisson's equation. This choice implies that 
variables and coefficients of the equation [T3] vary linearly in each element. For each node i 
of the system, the basis function is written 0j(-2). Then, the equation which has to be solved 
can be written as, 

^{z)dz = 0. 



+ 



d 



dz 



PwQL-ek I o 1 



(14) 



Using this method, all variables and coefficients are approximated on each element using 
the same basis function: 



1 

ne 

h{z,t) = ^hi{t)(l)i{z), 
1 

ne 

Cekiz,t) = y^^Cekiit)(pi{z), 

1 

ne 

ar{z,t) = ^ari{t)(j)i{z), 



(15) 
(16) 
(17) 
(18) 



10 V. Allegre, F. Lehmann, P. Ackerer, L. Jouniaux and P. Sailhac 

where nn = ne + 1 is the number of nodes in the system of length /, and ne the number 

of elements. The V^(t), hi{t), Cekiif) and ari{t) are respectively the values of the electrical 

potential, water pressure, electrokinetic coupling and bulk electrical conductivity at node i. 

For notations simplicity Ceki{t) will be written £j(t) in the following. After integrating by 

parts equation [141 the final system of equation to be solved is given by, 

[A].V, = -[B]X-{F,}, (19) 

where the elements of [A], [B] and vector {Fj} can be deduced from the following integrals: 

Aj = I 5^ 0fc CTfc J (P-^'jdz, (20) 

Bij = pw9 / f ^0fc £fc j (l)'i(p'jdz, (21) 

/I / nn \ 

ij^'t'kCkU'jdz. (22) 



'he matrices [A 



( jPress et al. 



is tridiagonal and the system can be solved using Thomas algorithm 



19921 ). The system introduced in equation [19] can be detailed for each element 



I as: 



[-{(Ti-i + (Ji).Vi-i + (cri_i + 2crj + ai+i).Vi 



2/\z 
-{ai + (Ti+i).V^i+i] = -^ [(/:i_i + 2Ci) + Ci+i).hi 

-{Ci-i + Ci).hi-i - {Ci + Ci+i).hi+i + {Ci-i + Ci+i)] (23) 

where A^ = 0.1 cm is the spatial discretization. Boundary conditions at nodes 1 and nn are: 

dV 
c^i^0i = -Jo (24) 

dV 

Crnn-K-<Pnn = -jl (25) 

with jo and ji are the values of the current density at the top and the bottom of the system 
respectively. The Poisson's equation is solved in terms of electrical potential V. Thus, two 
boundary conditions are necessary: one on the total current density J (Neumann type) and 
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one on the electrical potential itself V (Dirichlet type). As no current outflow can occur at 

the top and the bottom of the system (i.e no electrical exchange between the medium and 
the air), the total current density jo at z = and z = I are: jo = j« = A.m~^. In addition, 
a constant value of electrical potential Vq has been chosen at the top of the system. This 
constant can be a reference value, so that Vq = V has been chosen for simplicity. 
The modelling process is carried out using the following protocol: 1. Considering hydrody- 
namic boundary conditions, the equation [1] is first solved, so that water pressures and cor- 
responding water contents are computed at each node and each time step; 2. The Poisson's 
equation is then solved using equation [191 which allows to compute the electrical potential 
at each node of the system. 3. Thus, the needed potential differences can be deduced from 
the computed electrical potential field and compared to measurements. The next section 
presents a test performed considering a linear water pressure profile in saturated conditions, 
i.e simulating a Darcy's experiment. The following sections go into the unsaturated case and 
present computed SP differences and measurements in the case of a drainage experiment. 
Several hypotheses on SP coefficient dependence on water content will be considered and 
compared. 

3 SATURATED FLOW MODELLING 

The first step is to test the code accuracy by modelling a saturated flow through a column 
of sand (/ = 1.16 meter height), i.e a Darcy's experiment. The complete scheme described 



above was applied, so that equations [Tl and [TOJ were solved. The water electrical 
u,,, an d saturated SP coefficient Cgat corresponding to this test are those of 



conductivity 



AUegre et al. 



( 120101 ) and are reported in Table [TJ A linear and steady-state water pressure profile was 
applied to the medium (Figure [TtL). Thus, boundary conditions at the top and the bottom 
of the column are constant water pressure head and given by: ho{t) = 20 cm and hi{t) = 60 



12 V. Allegre, F. Lehmann, P. Ackerer, L. Jouniaux and P. Sailhac 

cm. In this case the medium remains saturated, so that S^^ = 1. The boundary conditions for 

Poisson's equation was a zero current density Jq at the top and a constant electrical potential 

Vo =0 V at the bottom of the column. The first simulated electrode (#10) is located at 11 

cm from the top of the system, while other electrodes are placed each 10 cm along it. The SP 

are then computed between the first five electrodes (#10 to #6) and the reference electrode 

(#1 15 cm above column's bottom) and for each dipole formed by two consecutive electrodes 

(e.g Al/io,9 will be the computed electrical potential difference between electrodes #10 and 

#9). The simulation were performed for 200 hrs, and the same results were obtained at each 

time step. 

The SP differences computed between each electrode and the reference are different. This is 

coherent with the increasing size of each dipole, so that larger the dipole is, larger the SP 

difference is (Figure [T|d). As the medium is homogeneous (fully saturated) SP differences 

computed for each dipole are merged (Figure [It) and exhibit the value AVij_i = —0.0954 

mV. The corresponding driving pore pressure difference for each dipole is APj i_i = 59.67 Pa. 

Considering these values, the inferred SP coefficient at saturation computed with equation 

[T^is: Csat = — 1.6x10^'^ V.Pa^^. This value is exactly equal to the original Csat (see Table 

[1]) used for calculation. This test experiment was performed using others water pressure 

profiles, other Csat values and other boundary condition value Vq and always yielded to 

the original Csat- We therefore confirm that the apparent measurement of AV/ AP for the 

dipoles can provide the streaming potential coefficient for saturated conditions. 
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4 DRAINAGE EXPERIMENT 

4.1 Unsaturated flow modelling 



We pr opose to model the SP measurements of a drainage experiment carried out by 



AUegre et al. 



( l2010l ). The SP were measured during a drainage experiment performed in a column of plexi- 
glass of approximately 1.2 metre height and 10 cm diameter fuUfilled with clean Fontainebleau 
sand. A constant water pr essure applied at the bottom of the column allowed the drainage to 



start. lAUegre et al. 



(120101 ) combined SP measurements to water content and water pressure 
measurements each 10 cm along the column. The aim of this section is to computed SP using 
the numerical scheme described above, considering sever al hypotheses for the S P coefficient 



depen dence on water content. Thus, formula for Cr of 



Guichet et al. 



(120031) 



Revil et al. 



(120071 ) and a new empirical model will be implemented to solve the Poisson's equation. 
At the beginning of the experiment, the medium is fully saturated, so that the water pres- 
sure profile is hydrostatic. The drainage starts when the boundary condition at the bottom 
of the column is set to a value oi h = 2 cm of equivalent water height (i.e Pw{z = /) ~ 
200 Pa). The water pressure heads and water contents are computed at each time step by 
solving the equation [T] and are compared to experimental data (Figure [2]). The parameters 
of the retention model (eq. [2]) and the relative permeability model (eq. [3]) are reported in 
Table [2l Since pressure sensors are located each ten centimetre along the column, before the 
drainage start, the pressure heads measurements are shifted from 10 cm between each other. 
These measurements are characteristics for the hydrostatic equilibrium. After the drainage 
start, the measured pressure heads decrease all at the same time and stabilize at different 
values depending on the sensor location. This shift of the pressure values at the end of the 
experiment (e.g at t~ 150 h) indicates the presence of a capillary fringe and of a water sat- 
uration gradient (see Figure |2]d). For very long time (around 90 days for the sand used) the 
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water phase equilibrates itself in the column until a linear water pressure profile is reached. 

From the drainage start, the water saturations do not decrease at the same time, but one 

after the other depending on the measurement location. The time shift between the water 

content measurements informs on the dynamic of the saturation front propagation during 

the drainage experiment. 

The hydrodynamic parameters used in the Richards' equation (Table [2]) were computed by 

inversion of the water pressure heads a nd water content m easurements. The detailed in- 



verse problem procedure can be found in 



Hayek et al. 



(120081 ). It is shown that this approach 



provides good estimations and small errors on 6'^, A and ha- 



4.2 Unsaturated SP modelling 



To solve the Poisson's equation at each time s t ep, a m odel for the SP coefficient dependence 



on water saturation is needed. 



Guichet et al. 



( I2OO3I ) proposed that the SP coefficient vary 



linearly with the effective water saturation. 



Or — Op_. 



(26) 



Revil et al 



(I2OO71 ) had a different approach and proposed another relation depending on a 



relative permeability model as. 



Or 



^-n+l 



With fcr = ^f+2+2/A 



(27) 



where n is the Archie's sat uration exponent (lArchie 



as L =0.5 (JMualem 



19421). The L parameter is usually chosen 



19761 ) , but is equal to 1 in iRevil et al 



( l2007l ). so that the two cases will 



be tested in the following. These two models have in common to predict the maximum of 
the SP coefficient to be Csat (for Sw = ^)- Moreover, they imply a mono tonous decrease of 



the relative SP coefficient with decreasing saturation. 



Allegre et al 



( I2OIOI ) recently observed 



that Cr could be until 200 times larger than Csat-, and that it first increases with decreasing 
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saturation, and then decreases up to the residual water saturation. So that, we propose an 

empirical relation for the SP coefficient inferred from these measurements written as. 



C = CsatSe[l + (3(1-8, 



(28) 



where f5 and 7 are two fitted parameters. Note that (3 depends on the considered dipole and 
varies as a function of the vertical location. This assumption coming from the experimental 
SP coefficients, which exhibit different maximum values, will be discussed in the following 
section. Contrary to the two ffist relations, this model predicts a non-monotonous behaviour 
of the SP coefficient as a function of water saturation. The three presented models were used 
to compute C{Se) from computed water saturations, to solve the Poisson's equation. 
Moreover, an a priori is needed for the electrical conductivity crr{Sw) to solv e equation [TOl 



The electrical conductivity is inferred at each time step using the Archie's law (jArchie 



19421), 



ar = aw(j)'^S^ 



(29) 



with m and n the two Archie's exponent, a^ the water electrical conductivity [S.m^^] and </> 
the porosity. The relation C{Se) is also computed at each time step and for each element of 
the system from computed water saturations. All parameters' values needed for computation 
are reported in Tables [2] and [31 Using the computed electrical potential values in the column, 
SP were computed for the dipoles (10,9), (9,8), (8,7) and (7,6), which are the dipoles located 
in the unsaturated part of the column during the drainage experiment. 



5 DISCUSSION 



The SP were computed s olving equations [J an d [TOJ 



The computed SP using 



Guichet et al 



f!2003h and 



or the three prese nted models (Figure |3]). 



Revil et al. 



( 120071 ) models are very small 



compared to the measurements (Figure [3tL,b). A jump, corresponding to the drainage start, 
is observed in SP signals at the beginning of the simulation. Its magnitude is around 0.09 
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mV for tests perfor med using; equat i ons [26] and [27] and is tlie same for all dipoles (Figure 



[3li,e). In the case of 



Guichet et al 



(120031) model (Figure [Sjl), the computed SP begin to 



decrease to reach a minimum value for dipoles (10,9) and (9,8) around -0.11 mV and then 



increase. This minimum is ob served at t 
for all dipoles. In the case of [Revil et al. 



l e tim e when water saturation stops to decrease 



( I2OO7I ) model, the increasing of computed SP is 



monotonous during the simulation, from a value around -0.09 mV to almost for all dipoles. 
It is obvious that these two models are not appropriate to explain the SP measurements both 
in terms of magni tude and behaviour. 



Revil et al. 



teOOTl ) used L =1 in their model, instead of L =0.5 as proposed by 



Mualem 



(ll976[ ) to hold as the best value for forty-five soils including sand. Cons equently, th e 



Revil et al. 



equa - 



feoOTh 



tion [27] was also implemented using L =0.5, which leads to a modified 

model. The results (Fig. [3^) are similar to those for L =1 in terms of amplitude, but show 

an increasing of SP signals without any step as it was observed before. The choice of L is 



quite important because it's invol yed in the glo 



influences the shape of the model ( lAllegre et al. 



3al p ower law in eq. [27] and consequently 



2010). 



On the other hand, the model introduced for C,- by equation [28] leads to good results in terms 
of computed SP. The measured SP signals are well reproduced, particularly at the begin- 
ning of the experiment between t ~ 20 h and t ~ 100 h. The computed SP corresponding to 
dipoles (10,9) to (7,6) decrease one after the other when water saturation (measured at the 
same level) begins to decrease (see Fig. [2}o and[3]3). Thus, the saturation front propagation 
is characterized by the time shift between SP decreasing starts. The computed SP using 
both equations [26] and [27] are up to one and a half order of magnitude smaller than those 
computed using expression [28] (Fig. [3]). 

For times over 100 hours, the residuals between measured and computed SP are larger. At 
this point of the experiment the water flow is very low, leading to very small variations of 
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SP. This point is interpreted in lAUegre et al. 



mdence on water content in sand 17 



(120101 ). who precise that for such low water 



outflow (even not measurable), the SP variations are too weak to insure robust interpreta- 
tion. Consequently, a better fit was expected at the beginning of the experiment (when the 
water fiow is maximum), which is the case for three measurement dipoles. In addition, all 
the tests performed to ensure the quality of SP recordings and prevent the external sources 
of noise from di sturbing the measur ements, including a statistical study on uncertainties, 



can be found in 



AUegre et al. 



(120101 ) (Appendices A, B). Moreover, the fit and estimation 
of the set of parameters of eq. [281 would be improved if the SP measurements were inverted. 
Thus a joint inversion would give more i nformations on para meters sensitivities. 



One can verify the assumption made by 



AUegre et al. 



(I2010h to use eq. (112]) to inferred SP 



coefficients from SP measurements. Thus, SP coefficients was computed using eq. (1121) . and 
compared to SP coefficients computed with eq. (128|) averaged on 10 cm to be representative 
of the same investigated volume (Fig. H]). It is shown that SP coefficients using eqs (IT2l) and 
( 128|) give very close results in terms of behaviour and magnitude. This example prove the 
validity of using eq. ( 1T2l) to infer accurate values of C even for non-steady conditions. One 
can consider that these values are apparent SP coefficients because the water saturation 
distribution is not perfectly homogeneous in 10 cm of sand, however they are still perfectly 
representative of true SP coefficients, at least for this kind of drainage experiment. 
Furthermore, computed SP coefficients using different assump tions for Cr (equatio ns l26l 1271 



AUegre et al. 



mm (Fig. ED . 



and 12S]) can be compared to measurements of AV/ AP from 
The fitted values of (3, 7, and parameters of equations [26] and [27] are reported in the Table 
[3] It is clear on figure [5] that the two existing models for Cr predict lower values of Cr than 
measurements and fail to provide the correct behaviour of C^. Moreover, the measured SP co- 
efficients AV/ AP are quite well re produced in the case of using equation [2H] It is important 



to notice that lAUegre et al.l (120 101 ) inferred experimental SP coefficients using the equation 
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[T2] which is classically used for steady-state flows, neglecting electrokinetic sources com- 
ing from SP coefficient and electrical conductivity contrasts. Since no electrokinetic sources 
have been neglected in the present modelling, one can conclude that the non-monotonous 
behaviour combined to large values of C,. observed previously are not artefacts and are not 
created by neglected contrasts of Cr or Oj. but has a physical origin. Therefore we show that 
the measurement of apparent AV/AP for the dipoles can provide the streaming potential 
coefficient Cr in these conditions. Nevertheless, a slight difference is observed between the 
maximum of computed and measured SP coefficients. These differences come from the use 
of the equation [12] to infer SP coefficients. Finally, these results confirm that a different 
maximum in C^ for each dipole is necessary to provide accurate computed SP differences. 
We suggest that this observation has a physical meaning coming from the flow behaviour. 
Indeed, each dipole is not affected by the same hydrodynamic conditions, and particularly 
the same flow velocity, when water saturation (at its level) decreases. 

An important point to discuss is the behaviour of the total current density. The flgure[6]shows 
six snapshots of both water saturation and current density component Jcond = —<Jr'^V and 
Jconv = cTrCVP Vertical proflles at six time during the simulation. The components Jcond 
and Jconv are computed a posteriori using flnite differences. Even if SP sources are linked to 
V- J, as introduced by the conservation equation (eq. [9]), we think that the vertical compo- 
nent of J is still relevant to describe the general behaviour of SP signals. This representation 
is usefull a posteriori to interpret the global SP variations, as it integrates all electrokinetic 
sources coming from SP coefficient or electrical conductivity contrasts (see eq. 10). The 
water saturation vertical profiles characterizes the propagation of saturation front during 
the drainage. In the saturated part of the column, 9 = Og thus S'u; = 1. In the unsaturated 
part, the saturation decreases from the top to the saturation front. The saturation value 
at the top of the system decreases during drainage from Sw{z = 0) = 0.5 at t = 6 h to 
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Sy^{z = 0) = 0.33 at t = 50 h and then stabilizes. For times greater than t = 50 h, the 

water saturation is almost constant (5*^ ^ 0.33) for < z < 40 cm. Moreover, it is shown 
on figure [6b that the conduction (Jcond) and the convection (Jconv) component of the total 
current density J are almost equal in absolute terms in the whole system during the entire 
experiment. The maximum difference observed between computed Jcond and Jconv does not 
exceed 2 per cent. This is a crucial point because the relation Jcond = Jconv allows the use 
of equation [12] (i.e C = AV/AP) to deduce SP coefficient values from measurements. This 
is another confirmation that this approach can be used and gives accurate results for this 
kind of experimental conditions. 

A large discontinuity is observed at the bottom of the column. This comes from the im- 
portant contrast of electrical conductivity between the saturated medium and the water in 
the reservoir at the bottom of the column. Thus, the electrical conductivity changes from 
asat — 0.002 S.m^^ to a^ = 0.01 S.m^^ at this boundary. This is the only possible source 
of current since the SP coefficient is zero in water. Nevertheless, this contrast does not in- 
fiuence the SP measurements, since the measurement are located far from it and that the 
current dens i ty ret urns to a constant value in the centimetre above it. It was suggested by 



Linde et al. 



(120071 ) and Revil and Linde {comment submitted), that this contrast could be 
responsible for the behaviour of the measured SP. However, it is not the case for our exper- 
iment but may occur for measurements involving an electrode very close to this source. It 
should be then considered as an experimental artefact in this case. 

The second point is that the maxima of Jcond and Jconv are not located at the higher gra- 
dient of water saturation (at the saturation front), corresponding to the higher contrast of 
electrical conductivity a^, but backward from this front. Its absolute value varies not linearly 
from 3x 10"^ A.m"^ to 6.8x 10"^ A.m~2 for 6 h < t < 117 h (Figure [7]). This suggests 
that the predominant contribution to the total electrical density comes from contrasts in 
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SP coefficient which is larger for S*^ ^ 0.8 than for saturated conditions. This is obviously a 

consequence of the behaviour of Cr as a function of water saturation. Thus, it seems that the 
electrokinetic response could be dominated by SP coefficient gradients and not by electrical 
conductivity gradients for such water flow. 

This statement can be investigated using different parameters in the Archie's law (eq. [29]) 
to increase the influence of the electrical conductivity. Thus, the previous value of n = 1.45 
is replaced by a larger value as n = 2.5. The influence of the saturation exponent n depends 
on the considered model for the SP coefficient (Fig. [8tL,b,c). Indeed, any change in n does 
not influence the computed SP in the case of using equations [26] or [28] for calculation. On 
the contrary, the increasing of n from 1.45 to 2.5 changes both the behaviour and magnitude 
of computed SP when the equation [27] is implemented (Fig. [Hb)- Thus, resulting SP exhibit 
larger values than the previous ones. This is explained by the increasing of the SP coefficient 
value implied by the increasing of n (see eq. [27]). It suggests that electrical conductivity 
contrasts are insignificant on the computed SP compared to other electrokinetic sources in- 
duced by SP coefficient contrasts. 

One can conclude that a non-monotonous behaviour combined to large values of C{Sw) needs 
to be implemented in the Poisson's equation to obtain computed SP close from the measured 
one. Nevertheless, the expression presented in this work is not sufficient and our approach 
does not deal with physical considerations. In future work, the water-flow conditions of the 
drainage experiment and especially the pressure dynamic and flow velocity could be involved 
in the observed SP coefficient behaviour. 

6 CONCLUSIONS 

We presented in this work the ffist modelling of both Richards and Poisson's equations for 
unsaturated conditions using ID finite element method. Several simulations have been per- 
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formed using three different hypotheses on the SP coefficient C r to de duce SP differences 



which have been compared to o bserv ations from 



models from IGuichet et aL 



(120031) and 



Revil et al 



AUegre et al 



(l2010l ). The two existing 



f|2007l) w e re not able to predict SP dif- 



AUegre et al 



(l2010l ). although no electrical 



ferences consistent with the measurements of 
current sources have been neglected. We proved with this modelling approach that a non- 
monotonous expression for C{Se) is necessary to correctly reproduced these measurements. 
We also demonstrated the consistency of our SP coefficient dataset verifying the equality 
Jcond = Jconv, and Consequently the possibility to use the apparent AV / AP measurements 
to infer correct C values for this kind of drainage experiment. This conclusion has to be 
verified for other flow conditions, for example steady-state unsaturated flow conditions with 
various flow velocities using different sands. 

Finally, a joint inversion approach of hydrodynamics and electrical potential, which could 
improve the modelling results of SP, is undergoing. This approach will help to insure the 
robustness of our model for long times and will define precisely the sensitivities of inverted 
model parameters. 



For field applications, even i 



and precipitation occurrence (JThony et al. 



a good correlation can be observed between SP response 



1997 



Doussan et al. 



2OO2I ). it seems that a lin- 



ear relationship between water flux and electrical potential gradient would be difficult to 
establish. Indeed, non-linear effects coming from the behaviour of the SP coefficient for un- 
saturated conditions, show that this relationship is probably more complex. Moreover, some 
acquisition i ssues and changing s oil conditions during long experiments make SP difficult 



to interpret fJDoussan et al. 



2OO2I ). However, shorter artificial infiltration experiments using 
SP monitoring could still be possible at the field scale. SP could be modelled with our new 
model of relative SP coefficient taking into account for infiltration and/or evaporation with 



time-varying upper boundary conditions. These experiments could be very useful to infer 
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some hydrodynamic parameters of soils, such as hydraulic conductivity, and give a way to 

measure ground water flux in the vadose zone. 
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Table 1. Parame ters used to perform the Darcy's experiment simulation, deduced from 



Allegre et al. 



(|20im . 



au, [S.m-i] Csat [V.Pa-i] ho [cm H2O] hi [cm H2O] 



103.2x10 



-4 



-1.6x10" 



20 



60 



Table 2. Hydrodynamic parameters values used to solve the Richards equation for unsaturated 
conditions. The parameter K^^"'^ is the measured permeability of the sand. These two values of 
permeability lead to the same results for pressure and water-content behaviours. 



Ks (xlO-5) [m.s-i] K^'"'' (xlO-5) [m.s^^] ha [m] A 6*,. [- 



bnf 



1.65 



17.2 



0.4 3.88 0.11 0.36 0.305 



Table 3. Parameters needed to implement equations [27l [28] and [29] in the Poisson's equation. The 
fou r values giveii for corresponds to the dipoles from (10,9) to (7,6). The n value was measured 



in (| Allegre et al. 



2OIOI). 



Csat [V.Pa ^] (Tw [S.m ^] n (Archie exponent) 



/3 



7 



-1.6x10^^ 103.2x10"^ 



1.45 



32,52,72,92 0.4 
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Figure 1. a) The water pressure profile used for the Darcy's experiment simulation. The pressure 
head boundary conditions are constant and equal to /iq = 20 cm and hi = 60 cm of equivalent 
water height at the top and the bottom of the column respectively, b) Computed SP between the 
first five electrodes (#10 to #6) and the reference #1. c) Computed SP for the dipoles (10,9), 
(9,8), (8,7) and (7,6). The computed values are AVi,i - 1 = -0.0954 mV and APi,i- 1 = 59.67 
Pa which lead to Csat = LGxlQ-^ V.Pa^i. 
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Figure 2. a) Measured (dots) and computed (dashed lines) water pressure iieads deduced from tlie 
Richards equation solving. Indices i indicates the location of measurement, as /iio is the pressure 
head 16 cm from the column top and hi 10 cm from the column bottom. Other h values are 
measured and computed each 10 cm. b) Measured (dots) and computed (dashed lines) water 
saturations deduced from the Richards equation solving, where indices i indicates the location of 
measurement. The drainage starts at t ~ 22 hr. 
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Figure 3. Measured (dots) and computed (lines) streaming potentials deduced from the Poisson's 
equation solving, for each dipole, using respectively equations [261 (a) . [271 (b) and [281(c). Computed 
SP using equ ations 1261 fd) . 1271 fe ) and[28](f) in the Poisson's equation, and parameters from tables 



[2] and m The 

(dashed lines). 
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Figure 4. Computed relative SP coefficients using eq. (J12p with computed AV and AP after 
Poisson's equation was solved (lines), and using eq. (j28p (dashed lines), for five locations in the 
column. 
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Figure 6. Snapshots of the water saturation profile (a) and total current density components Jcond 
(black line) and Jconv (I'ed line, b) for six times between t = 6 h and t = 117 h. The electrical 
current density is expressed in A.m~^. The snapshot corresponds to the simulation performed using 
the equation [28] for implementation of C{Se) in the Poisson's equation. 
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Figure 7. Maximum of the total current density J (black line) inferred from curves in Figure [6] 
and corresponding water flow velocity (dashed black line). 
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Figure 8. Computed SP (lines) using equations [26] (a) . [271 f b) and [281(c) with n = 1.45, compared 
to computed SP with n = 2.5 (dashed lines). 
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